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A  new  dynamical  core  for  numerical  weather  prediction  (NWP)  based 
on  the  spectral  element  Eulerian-Lagrangian  (SEEL)  method  is  presented. 
This  paper  represents  a  departure  from  previously  published  work  on  solv¬ 
ing  the  atmospheric  equations  in  that  the  horizontal  operators  are  all  writ¬ 
ten,  discretized,  and  solved  in  3D  Cartesian  space.  The  advantages  of 
this  new  methodology  are:  the  pole  singularity  which  plagues  all  grid- 
point  methods  disappears,  the  horizontal  operators  can  be  approximated 
by  local  high-order  elements,  the  Eulerian-Lagrangian  formulation  permits 
extremely  large  time-steps,  and  the  fully-implicit  Eulerian-Lagrangian  for¬ 
mulation  only  requires  the  inversion  of  a  2D  Helmholtz  operator.  In  order  to 
validate  the  SEELAM  model,  results  for  four  test  cases  are  shown.  These 
are:  the  Rossby-Haurwitz  waves  number  1  and  4,  and  the  Jablonowski- 
Williamson  balanced  initial  state  and  baroclinic  instability  tests.  Compar¬ 
isons  with  four  well-established  operational  models  show  that  SEELAM  is 
as  accurate  as  spectral  transform  models. 


1.  INTRODUCTION 

Because  of  the  changing  trends  in  high  performance  computers  from  large  vector  machines  to 
distributed-memory  architectures,  numerical  methods  that  decompose  the  physical  domain  into 
smaller  pieces  have  been  receiving  significant  attention.  This  new  focus  on  local  methods  is  es¬ 
pecially  true  in  the  atmospheric  sciences  where  very  large  models  covering  the  entire  globe  are  run 
in  time-scales  ranging  from  days  (in  numerical  weather  prediction)  to  thousands  of  years  (in  cli¬ 
mate  simulations).  Finite  difference  and  finite  element  methods  are  two  such  methods  which  decom¬ 
pose  the  domain  locally  thereby  facilitating  their  implementation  on  distributed-memory  computers. 
However,  one  of  the  biggest  disadvantages  of  these  methods  is  that  traditionally  they  have  not  been 
able  to  compete,  in  terms  of  accuracy,  with  spectral  transform  methods  which  are  typically  used 
operationally  in  numerical  weather  prediction  (NWP)  and  climate  modeling.  For  example,  spectral 
transform  models  are  used  by  the  National  Center  for  Environmental  Prediction  [8],  the  European 
Centre  for  Medium-Range  Forecasts  [9],  the  National  Center  for  Atmospheric  Research  [3],  and  the 
U.S.  Navy  [4]. 

Spectral  element  methods  combine  the  local  domain  decomposition  property  of  finite  element 
methods  with  the  high-order  accuracy  of  spectral  transform  methods.  In  other  words,  spectral 
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elements  are  as  local  as  finite  element  methods,  and  thereby  can  be  used  as  efficiently  on  distributed- 
memory  computers  while  sustaining  the  same  level  of  accuracy  obtained  with  spectral  transform 
methods.  These  methods  are  essentially  high-order  finite  element  methods  where  the  grid  points 
are  chosen  to  be  the  Legendre-Gauss-Lobatto  (LGL)  points.  This  choice  of  grid  points  allows  for 
stable  high-order  interpolations  and  results  in  efficient  numerical  integration  strategies  because  the 
LGL  points  are  also  used  as  the  quadrature  points  in  the  numerical  integration  required  by  the  weak 
integral  formulation  common  to  all  Galerkin  methods. 

Using  the  spectral  element  method  we  have  developed  a  hydrostatic  primitive  equation  model 
for  NWP  [2].  However,  in  that  work  the  time  integration  scheme  used  was  an  explicit  Eulerian 
method.  Although  explicit  time-integration  methods  for  atmospheric  simulations  can  be  incredibly 
fast  and  quite  accurate,  their  main  disadvantage  is  that  small  time  steps  must  be  observed  in  order 
to  maintain  stability.  The  reason  for  this  prohibitively  small  time-step  is  due  to  the  fast  moving 
gravity  waves.  These  waves  require  a  small  time-step  while  only  carrying  a  very  small  percent  of 
the  total  energy  in  the  system.  In  order  to  ameliorate  this  rather  stringent  time-step  restriction 
atmospheric  scientists  have  tried  various  approaches  such  as  using  a  larger  differencing  stencil  for 
the  gravity  wave  terms  thereby  effectively  reducing  the  Courant  number,  and  discretizing  the  gravity 
wave  terms  implicitly  in  time.  We  can  easily  employ  the  first  strategy  in  our  current  formulation, 
that  is  using  a  larger  differencing  stencil  for  the  gravity  wave  terms,  but  this  is  more  typically 
done  to  avoid  the  inf-sup  (Babuska-Brezzi)  condition.  However,  discretizing  the  gravity  wave  terms 
implicitly  in  time  is  a  much  more  effective  way  of  increasing  the  time-step. 

After  the  gravity  wave  terms  have  been  successfully  discretized  the  next  set  of  terms  responsible  for 
controlling  the  maximum  time-step  are  the  advection  terms.  In  order  to  use  increasingly  larger  time 
steps  atmospheric  scientist  have  turned  to  Lagrangian  methods  for  treating  these  recalcitrant  terms. 
By  rewriting  the  equations  in  Lagrangian  form  the  troublesome  advection  terms  are  absorbed  into 
the  total  derivative.  Thus  the  equations  in  this  form  are  now  discretized  in  time  along  characteristics 
which  results  in  a  much  more  stable  numerical  method  due  to  the  disappearance  of  the  Courant- 
Friedrichs-Lewy  (CFL)  condition. 

In  the  current  paper  we  combine  a  Lagrangian  method  with  high  order  basis  functions  as  we  showed 
in  [1]  and  extend  this  hybrid  method  to  the  solution  of  the  hydrostatic  atmospheric  equations.  The 
allure  of  this  method  is  that  it  achieves  the  same  order  of  accuracy  obtained  with  exponentially  high 
order  explicit  Eulerian  methods  [2]  while  permitting  time-steps  at  least  5  times  larger. 

2.  ATMOSPHERIC  EQUATIONS 

In  this  paper  we  show  how  to  solve  the  hydrostatic  primitive  equations  which  describe  the  motion 
of  the  atmosphere.  We  assume  dry  physics  (i.e.,  no  physical  parameterization)  and  thus  only  take 
into  account  the  dynamical  processes.  The  equations  we  solve  are 


du 

dt 


dn  8 

w  +  v.(™)  +  -(^)  =  0 


.  du  2 ilz  8P 

U  ■  VM  +  a— —  = - X  U)  —  \  Lp  —  C„P— —  V 7T 

da  a -  ott 


89  ^  .89  n 

—  +u- V9  +  a—  =  0 
dt  da 


8P  =  -Cp° 
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where 


7T  =  Ps  -  Pt  ,  &  = 


P~Pt 


7 r 


and  (a.  Cl)  are  the  radius  and  angular  rotation  of  the  earth,  respectively.  The  prognostic  variables 
for  these  equations  are  q  =  (ir1u,9)T,  while  a  and  ip  represent  diagnostic  variables. 

3.  TEMPORAL  DISCRETIZATION 

The  Lagrangian  form  of  the  conservation  equations  are 


dir  (_  da 

r-T,+s; 


du  2Clz  .  „dP _ 

—  = - o~{x  xu)-V<p-  Cp8—-Vir 

at  az  air 


?=0 

dt 


(5) 

(6) 
(7) 


where  the  Lagrangian  derivative  is 


d  d  „ 


Discretizing  the  conservation  equations  in  time  by  the  second  order  backward  difference  formula 
yields 

i^\  n+1 

(8) 


7T"+1  -  §7Tn  +  ITT”"1 


At 


2  (ri  d°X 
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un+1  -  + iw"-1 

At 


2  Clz  f)P 

— 2~(x  X  u)  +  +  CpO^—V 

az  077 


n+1 
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At 


where  the  variables  qn  and  qn  1  denote  Lagrangian  departure  point  values. 


(9) 

(10) 


3.1.  Linearization  of  the  Equations 

To  avoid  having  to  solve  a  nonlinear  problem  we  linearize  the  equations  as  follows.  The  surface 
pressure  is  linearized  about  the  reference  state  7r*  =  1000  hPa  as  such 


7T  =  7 T  —  7T 


where  7r  now  denotes  the  perturbation  surface  pressure.  In  addition,  let  us  introduce  the  reference 
potential  temperature 

rp* 

A*  —  1 

~  P(tt*) 

with  T*  =  300°  Kelvin.  It  should  be  understood  that  the  reference  potential  temperature  is  a 
function  of  the  vertical  coordinate  er  because  the  Exner  function  varies  with  cr. 


3.1.1.  Surface  Pressure  Equation 
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With  the  above  linearizations  the  surface  pressure  equation  takes  the  following  form 


7r"+1  -  f  TT71 


l-n-i 
3  H 


At 


— - 7 r 

3 


_  da\ 

v-u  +  ai) 


n+1 


(ii) 


3.1.2.  Hydrostatic  Equation 

Beginning  with  the  finite  differenced  equation  for  geopotential 


‘fik  —  (Pk+1  —  Cp9k  (Pk+ 1/2  ~  Pk )  +  Cp9k+ 1  (Pk+ 1  —  -Cc+l/2) 

where  k  =  l,...,Niev  and  TVjet,  represents  the  number  of  vertical  levels  in  the  model,  we  can  then 
take  a  Taylor  series  expansion  about  the  reference  surface  pressure  7r*  and  potential  temperature  6* 
yielding 


<Pk  fk+ 1  —  CpOk  {Pk+1/2  Pk )  +  cp@k+l  (Pk+1  Pk+l/2 


(12) 


cv0l  \ 


dPk+l/2  9P£\  ,  ,  n*  (dPk+ 1  dPk+ 1/2  \, 

1  (7T  -  7T  )  +  CpVk+1  (  — ^ - 7^ -  I  (7T  -  7T  )  . 


dn  )  "  '  '  ^“K+1  V  97T 

Equation  (12)  can  be  written  in  the  matrix  form 

Ak,m  =bk(e)  +  ck* 

where 

bm  =  cp0k  (p*k+1/2  -p;)+  CpOk+1  (pk+i  -  Pfe*+1/2) 


and 


—  r  a* 

Ck  ~  CP°k 


dPk+ 1/2  dP,*\  .  ...  ...  (dP\ 


(7T  -  7T*)  +  Cp9*k+1 


di r  dir 

This  results  in  the  following  geopotential  gradient 


r)P* 

k+ 1  UJrk+ 1/2 


<97T 


97T 


(7T  ^  7T*)  . 


(13) 


5.1.5.  Momentum  Equations 

Linearizing  the  gradient  due  to  surface  pressure  yields 


„  aPfc  A"+1  (  n*dPk^n+l 

cP9k  9jt  Vtt)  -  (cp9k  ^  Vtt 


The  momentum  equations  then  look  as  follows 


n+1  2  (2  Viz, 

u  + 1  +  -At  — —(a;  x  it) 
3  \  az 

which  in  matrix  form  are 
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Substituting  the  linearized  geopotential  from  Eq.  (13)  gives 


/)«».«+!  _  +  .n  _  l.n- 1 
A  Uk  ~  £Uk  ,juk 


'-At 


v^)-1  [bf(en+1)  +  cfnn+1])+cpe, 


r)P* 

*  Urk  y^r^n+1 

dn 


and  inverting  yields 


un+l  = 


((^iI)-i&r(en+i)) 

V  ({All)~1cf*n+1)  +cp^^V7f”+1 


M  -  K""1  - 


Let  us  now  rewrite  this  more  compactly  as  follows 


u 


n+1  =  (A™)-1  [buk  (r+1)  -  |a tel  (tt"+1)) 


where 


and 


K  (en+1)  =  Ul  -  -  ?At  [v  ((^;)-1&r(^+1) 


c£  (7T-+1)  =  V  ((^,«)-1cr  t"+1)  +  cP9l^Vn 


n+l 


(14) 


(15) 


(16) 


3.2.  Calculation  of  Departure  Point  Values 

The  discretization  of  the  hydrostatic  primitive  equations  that  we  have  described  so  far  are  complete 
provided  that  we  have  a  means  of  computing  the  state  vector  at  the  departure  point  values.  Below 
we  describe  two  possible  approaches  which  are  rather  similar:  the  semi-Lagrangian  method  and  the 
Eulerian-Lagrangian  method. 

3.2.1.  Semi-Lagrangian  Method 

In  the  semi-Lagrangian  method,  the  departure  point  values  are  obtained  as  follows.  We  solve  the 
Lagrangian  pure  advection  problem 


dq  dx 

—  =  0  and  =  u. 
dt  dt 


(17) 


In  the  current  implementation  we  use  a  2nd  order  Runge-Kutta  method  to  solve  the  trajectory 
equation. 


3.2.2.  Eulerian-Lagrangian  Method 

Although  this  method  is  known  in  the  literature  as  the  operator-integration-factor  splitting  method 
(see  [6])  we  have  chosen  to  refer  to  it  as  the  Eulerian-Lagrangian  method  because  it  represents 
an  Eulerian  approximation  to  the  Lagrangian  problem.  In  the  Eulerian-Lagrangian  method  the 
departure  point  values  qn  and  qn~ 1  are  obtained  as  follows.  We  solve  the  Eulerian  pure  advection 
problem 


Dqk 

Or 


+  u  ■  Vqk  =  0 


(18) 
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with  q(x,tn+1~k)  =  q(x,tn+1~k)  where  the  velocity  field  u  is  extrapolated  from  the  known  velocity 
fields  at  times  n  —  2,  n  —  1,  and  n.  In  the  current  implementation  we  use  a  2nd  order  Runge-Kutta 
method  to  solve  this  equation. 


4.  HELMHOLTZ  OPERATOR 

In  order  to  solve  for  the  surface  pressure  we  need  to  integrate  Eq.  (11)  vertically  from  a  =  0  to 
cr  =  l  which  yields 


JVleu 

E 


i 


At 


Aak 


N,ev 


E(v-<+1 


A  Ufe  - 


cr— 1 
cr— 0 


where  7r"+1  is  purposely  not  given  a  vertical  level  subscript  to  denote  that  the  arrival  surface  pressure 
is  only  a  function  of  the  2D  surface.  Using  the  no-flux  boundary  condition  at  the  top  and  bottom 
of  the  atmosphere  yields 


Nlev 

E 

fc=i 


n+1  4  Cin  ,  ljsra— 1 


7 r  1  —  o7 r, 


Me 


3;‘fc  T-  3  nk 


At 


-A ak  =  --i*  ^  (V  •  «)5!+1  Acrfc. 


(19) 


fc=i 


Substituting  un+1  from  Eq.  (14)  and  rearranging  gives 


NUv 


Me 


E*"+'a«  =  E 


1, 


-  -Atn*V 


k- 1  fc=l 

Substituting  for  cj(  (7?)  from  Eq.  (16)  yields 


(A”)-1  (  (en+1)  -  -Afcfeu  (tt”+1) 


Aon. 


Me„ 

E 

fc=i 

Me„ 


fc  =  l 


7fn+1 
4 


-  -AtW  ■ 


(A")-1  (v  ((A^)_1cf  7Tn+1)  +cp6»; 


dF\ 


dn 


k  V7f”+1 


Aak  =(20) 


E  l?- r^T1  -  xAfTr* V  •  ((Au)~1bk  (dn+1))  Aak. 


The  interesting  result  from  this  linearization  is  that  the  terms  (A^;)_1c^  and  cp0Jc  are  not 
functions  of  the  horizontal  operators  and  so  they  can  be  factored  from  the  gradient  and  divergence 
operators.  Doing  so  yields 


7f”+1  - 


Met, 


Me 


fc= 1 


-AfV  ( Ak,l )  V  +C; 


Puk  ' 


.dPl 

dn 


Act  k 


V  •  ((A“)-1V 


-lY7^«+l) 


(21) 


Y,\rk~  ^r1  -  lAt**v  ■  {(AurlK  (^n+1)) }  a^ 

k= 1  ^  ' 

where  we  have  removed  7rn+1  from  the  summation  because  it  is  not  a  function  of  the  vertical 
coordinate  <7.  This  expression  tells  us  that  the  first  bracketed  term  on  the  LHS  can  be  computed 
once,  because  it  is  not  a  function  of  time,  and  then  taken  as  a  constant  in  the  solution  of  the  2D 
surface  pressure  Helmholtz  equation.  To  this  end,  we  rewrite  the  Helmholtz  equation  as  follows 


7fn+1  -  AV  •  ((A“)-1V7rn+1)  =  b 71 


(22) 
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where 

4  /  r)P*\ 

A = n- AiV  E  w1# + Afffc  (23) 

k- 1  '  ' 

and 

N/eu  f  4  1  _  O  'l 

b*  =  E  3^3^_1  -  3 V  •  ((^T  ^  (O)  A^-  (24) 

fc= l  ' 

5.  GRID  GENERATION  ON  THE  SPHERE 

Hexahedral  (a.k.a.  cubic  gnomonic)  grids  are  constructed  by  subdividing  the  6  faces  of  a  hexa¬ 
hedron  into  the  desired  number  of  quadrilateral  elements,  and  then  mapping  these  onto  the  sphere. 
This  approach  results  in  the  construction  of  a  hexahedral  grid  with  the  following  properties 

Np  =  6(nHN)2  +  2 

Ne  =  6  (nn)2 

where  Np  and  Ne  denote  the  number  of  points  and  elements  comprising  the  grid.  The  parameter 
nn  refers  to  the  number  of  quadrilateral  elements  in  each  direction  contained  in  each  of  the  6  initial 
faces  of  the  hexahedron  and  N  is  the  order  of  the  basis  functions  of  each  of  the  elements.  The 
comparable  hexahedral  resolution  to  the  spectral  triangular  truncation,  T,  can  be  obtained  by  the 
expression 

T  =  hhN. 

Thus  to  obtain  T160  we  can  use  tin  =  20  and  N  =  8. 

6.  RESULTS 

6.1.  Rossby-Haurwitz  Waves 

In  order  to  judge  the  accuracy  of  SEELAM  we  have  run  Rossby-Haurwitz  waves  numbers  one 
and  four  and  compared  the  surface  pressure  of  SEELAM  with  those  of  NOGAPS  [4]  for  a  T160  L24 
resolution  after  5  day  integrations.  The  SEELAM  model  uses  a  time  step  three  times  larger  than 
the  NOGAPS  model  for  these  tests. 

Figure  1  shows  the  wave  one  results  while  Fig.  2  shows  the  wave  four  results.  There  are  slight 
differences  in  the  shape  of  the  waves  but  both  models  yield  the  same  maximum  and  minimum  contour 
levels.  More  importantly  both  models  yield  the  same  phase  speeds. 

6.2.  Jablonowski-Williams  Tests 

The  following  two  cases  represent  a  new  set  of  tests  for  judging  the  accuracy  and  stability  of 
dynamical  cores.  These  tests  are  introduced  in  Jablonowski  and  Williamson  [5]. 

6.2.1.  Balanced  Initial  State 

For  this  test  case,  the  atmosphere  is  initially  balanced.  With  the  given  initial  conditions  the 
equations  should  remain  balanced  for  an  indefinite  amount  of  time.  Figure  3  shows  the  normalized 
surface  pressure,  n,  L2  error  norm  as  a  function  of  time  for  a  four  week  period  for  SEE- AM  and 
SEELAM  with  T160  horizontal  resolution  and  24  vertical  levels.  SEE-AM  is  the  explicit  Eulerian 
version  presented  in  [2].  Note  that  while  the  error  oscillates  with  time  it  remains  bounded  which 
confirms  that  the  initial  balanced  state  is  maintained  by  both  models.  In  addition,  both  models  give 
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FIG.  1.  Rossby-Haurwitz  Wave  Number  4:  The  surface  pressure  (hPa)  for  a)  NOGAPS  and  b)  SEELAM  for 
T160  and  ATiev  =  24  for  a  5  day  integration. 
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FIG.  2.  Rossby-Haurwitz  Wave  Number  4:  The  surface  pressure  (hPa)  for  a)  NOGAPS  and  b)  SEELAM  for 
T160  and  N\ev  =  24  for  a  5  day  integration. 


identical  errors  up  to  18  days  at  which  point  the  models  differ  but  not  significantly.  This  result  is 
perhaps  the  most  encouraging  result  in  the  suite  because  it  confirms  that  the  Lagrangian  implicit 
SEELAM  is  behaving  like  the  explicit  Eulerian  SEE- AM. 

6.2.2.  Baroclinic  Instability 

This  case  is  similar  to  the  balanced  initial  state  except  that  now  a  perturbation  is  added  to  the 
initial  zonal  velocity.  This  perturbation  grows  until  a  baroclinic  instability  develops  and  then  breaks 
near  day  nine.  Figure  4  shows  the  minimum  surface  pressure  ps  as  a  function  of  time  for  SEELAM 
against  various  models  including  the  NCAR  spectral  transform  model  [3],  the  NASA  Goddard  finite 
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FIG.  3.  Jablonowski- Williamson  Balanced  Initial  State:  The  normalized  surface  pressure,  7 r,  L2  error  as  a 
function  of  days  for  SEE- AM  (dashed  line)  and  SEELAM  (solid  line)  for  T160  (njj  =  20  and  N  =  8)  with  24  vertical 
levels. 


volume  model  [10],  and  the  German  Weather  Service  icosahedral  finite  difference  model  [7]  which 
we  denote  as  GME;  the  results  of  the  latter  three  models  are  courtesy  of  Christiane  Jablonowski. 
The  results  in  Fig.  4  are  summarized  as  follows.  SEELAM  compares  well  with  the  three  established 
models.  The  lower  order  NASA  and  GME  models  give  very  similar  results.  Finally,  the  higher  order 
NCAR  and  SEELAM  models  compare  extremely  well  with  each  other. 

7.  CONCLUSION 

A  new  dynamical  core  for  numerical  weather  prediction  (NWP)  based  on  the  spectral  element 
Eulerian-Lagrangian  (SEEL)  method  is  presented.  In  a  previous  paper  [2]  we  showed  the  advantages 
of  using  spectral  elements  in  3D  Cartesian  coordinates.  In  this  paper  we  have  extended  the  explicit 
Eulerian  method  described  in  [2]  to  an  implicit  Lagrangian  method.  The  advantage  of  using  the 
implicit  Eulerian-Lagrangian  method  is  that  it  permits  time  steps  5  times  larger  than  the  explicit 
Eulerian  method.  This  increase  in  permissible  time  step  should  translate  into  increased  efficiency 
of  the  model.  While  the  SEELAM  model  has  not  yet  been  fully  tested  it  has  passed  its  first  four 
tests  demonstrating  that  it  behaves  similarly  to  well-established  climate  and  weather  prediction 
models.  These  models  include:  the  U.S.  Navy’s  NWP  model  and  the  NCAR  climate  model.  In 
order  to  become  competitive  with  these  well-established  models  the  iterative  solvers  for  inverting 
the  resulting  Helmholtz  operator  must  be  fully  optimized.  This  is  the  topic  of  future  work. 
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FIG.  4.  Jablonowski- Williamson  Baroclinic  Instability:  The  minimum  surface  pressure  (hPa)  as  a  function 

of  days  for  the  NASA  (finite  volume),  GME  (finite-difference),  NCAR  (spectral  transform),  and  SEELAM  (spectral 
element  with  nn  =  20  and  N  =  8)  models  using  26  vertical  levels.  (The  data  for  the  first  three  models  are  courtesy 
of  Christiane  Jablonowski.) 


REFERENCES 

1.  F.X.  Giraldo  and  J.B.  Perot,  A  spectral  element  semi-Lagrangian  method  for  the  spherical  shallow  water  equations, 
Journal  of  Computational  Physics  in  review  (2003). 

2.  F.X.  Giraldo  and  T.E.  Rosmond,  A  scalable  spectral  element  Eulerian  atmospheric  model  (SEE-AM)  for  NWP: 
dynamical  core  tests,  Monthly  Weather  Review  in  review  (2003). 

3.  J.J.  Hack,  B.A.  Boville,  B.P.  Briegleb,  J.T.  Kiehl,  P.J.  Rasch,  and  D.L.  Williamson,  Description  of  the  NCAR 
community  climate  model  (COM2),  NCAR  Technical  Note  NCAR/TN-382+STR.  (1992). 

4.  T.F.  Hogan  and  T.E.  Rosmond,  The  description  of  the  Navy  Global  Operational  Prediction  System’s  spectral 
forecast  model,  Monthly  Weather  Review  119,  1786-1815  (1991). 

5.  C.  Jablonowski  and  D.L.  Williamson,  Baroclinic  instability  test  with  two  jets  in  the  midlatitudes,  Abstracts,  The 
2002  Workshop  on  the  Solution  of  Partial  Differential  Equations  on  the  Sphere,  Toronto,  Canada,  10-10  (2002). 

6.  Y.  Maday,  A.T.  Patera,  and  E.M.  Ronquist,  An  operator-integration-factor  splitting  method  for  time-dependent 
problems:  applications  to  incompressible  fluid  flow,  Journal  of  Scientific  Computing  5,  263-292  (1990). 

7.  D.  Majewski,  D.  Liermann,  P.  Prohl,  R.  Ritter,  M.  Buchhold,  T.  Hanisch,  G.  Paul,  W.  Wergen,  and  J.  Baumgard¬ 
ner,  The  operational  global  icosahedral-hexahedral  gridpoint  model  GME:  description  and  high-resolution  tests, 
Monthly  Weather  Review  130,  319-338  (2002). 

8.  J.G.  Sela  J.G.,  Spectral  modeling  at  the  National  Meteorlogical  Center,  Monthly  Weather  Review  108,  1279-1292 
(1980). 

9.  A.J.  Simmons,  D.M.  Burridge,  M.  Jarraud,  C.  Girard,  and  W.  Wergen,  The  ECMWF  medium-range  predic¬ 
tion  models  development  of  the  numerical  formulations  and  the  impact  of  increased  resolution,  Meteorology  and 
Atmospheric  Physics  40,  28-60  (1989). 

10.  K.S.  Yeh,  S.J.  Lin,  and  R.B.  Rood,  Applying  local  discretization  methods  in  the  NASA  finite-volume  general 
circulation  model,  Computing  in  Science  and  Engineering  4,  49-54  (2002). 


